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Abstract 


Steam generator tube rupture (SGTR) accident is an important scenario needed 
to be considered in the safety analysis of lead-based fast reactors. When the steam 
generator tube breaks close to the main pump, water vapor will enter the reactor core, 
resulting in a two-phase flow of heavy liquid metal and water vapor in fuel assemblies. 
The thermal-hydraulic problems caused by the SGTR accident may seriously threaten 
reactor core’s safety performance. In this paper, the open source CFD calculation 
software OpenFOAM was used to encapsulate the improved Euler method into the 
self-developed solver LBEsteamEulerFoam. By changing different heating boundary 
conditions and inlet coolant types, the two-phase flow in the fuel assembly with 
different inlet gas content was simulated under various accident conditions. The 
calculation results show that the water vapor may accumulate in edge and corner 
channels. With the increase of inlet water vapor content, outlet coolant velocity 
increases gradually. When the inlet water vapor content is more than 15%, the outlet 
coolant temperature rises sharply with strong temperature fluctuation. When the inlet 
water vapor content is in the range of 5% to 20%, the upper part of the fuel assembly 
will gradually accumulate to form large bubbles. Compared with the VOF method, 
Euler method has higher computational efficiency. However, Euler method may cause 
an underestimation of the void fraction, so it still needs to be calibrated with future 
experimental data of the two phase flow in fuel assembly. 
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1 Introduction 


In 2002, six advanced nuclear energy systems including lead-based fast reactors 
(LFRs) were formally established as Gen-IV reactors!'!, The structure of lead-based 
fast reactors has been significantly simplified compared to that of other reactors!!. 
Lead-based fast reactors have advantages, such as good nuclear fuel transmutation 
capability, excellent economics, and inherent safety?!. Therefore, many countries have 
invested significant research and development efforts. China approved the China 
Initiative Accelerator Driven System (CiADS) in 2015, The CiADS consists of three 


5-6]. The reactor is a subcritical 


subsystems: accelerator, spallation target, and reactor! 
fast neutron reactor cooled by the liquid lead—bismuth alloy. The steam generator (SG) 
is the kernel component of the core heat transfer in reactor systems. LFRs usually 
adopt a pool structure design in which the core, main pumps, and steam generators are 
immersed together in the primary vessel coolant!. This design separates the 
pressurized water in the secondary loop from the liquid lead—bismuth only through the 
wall of the heat transfer tube of the steam generator. Simultaneously, many tube 
bundles were implemented to enhance the heat exchange capability inside the tube 
bundles. Moreover, there are large pressure and temperature differences on both sides 
of the tubes, which may cause significant thermal and mechanical stresses on the tube 
wall. Owing to the vibration and corrosion of the coolant in the first and second 
circuits, the heat transfer tube may be the weakest part of the primary circuit system". 
Therefore, the possibility of steam generator tube rupture cannot be ignored. 

When the steam generator pipe breaks, the high-pressure subcooled water in the 
secondary circuit is injected into the low-pressure and high-temperature liquid lead 
bismuth in the primary circuit. High-pressure subcooled water quickly vaporizes into 
steam, thereby forming a two-phase liquid metal-steam two-phase flow]. Currently, 
international research on steam generator tube rupture (SGTR) accidents in lead-based 
reactors is at an early stage. Relevant experimental studies have been conducted on 
the interaction mechanism between water and liquid metals. Sibamoto et al. 


designed a probe that can simultaneously measure the temperature and phase state of 


a liquid metal—water—vapor multiphase flow. The cavity formation process of water in 
liquid metal, cavity surface mode boiling, and temperature distribution of the 
multiphase flow were investigated. Sa et al.!'!! built an experimental device for 
dripping molten liquid metal into water and studied the steam explosion reaction 
mechanism between water and a lead—bismuth alloy in an SGTR accident. Later, Sa et 
al.l] used an alcohol solution to simulate a water medium and a fluorinated liquid 
working medium to simulate liquid lead—bismuth and conducted an experiment on 
alcohol injection into a high-temperature fluorinated liquid. They also simulated the 
boiling phenomenon of direct contact between water and high-temperature liquid 
metal in an SGTR accident. Dostal et al.!'?! built a direct-contact boiling circuit for 
water and liquid lead—bismuth and studied the heat transfer behavior of two-phase 


1415] conducted experiments on the contact fragmentation of 


boiling. Huang et al.! 
high-temperature lead—bismuth droplets and liquid columns with subcooled water on 
a platform where molten metal and water interact. Deng et al.!'®! conducted numerous 
experiments to inject water lumps into a molten lead pool at the Sun Yat-sen 
University. 

Owing to the radiation-shielding property of lead—bismuth, the image of the 


two-phase flow pattern detected by gamma rays is unclear!!”!, Therefore, numerical 


calculations were used to simulate two-phase and multiphase flow problems in SGTR 
accidents. Wang et al.l!3] used the safety analysis program SIMMER-IIl to evaluate 


the pressure evolution of a steam generator tube rupture accident, the migration of 


water vapor in the primary circuit, and the possibility of entering the core. 
Ciampichetti et al.l?! used SIMMER-IIl to simulate the interaction between 


high-pressure water and liquid lead during an SGTR accident and conducted a 
preliminary assessment of the pressure change of the covering gas in the vessel and 
the possibility of water vapor entering the core after the accident. Zhixing et al.!?°! 
simulated a breakdown accident of the upper and lower heads of a heat exchanger in a 
small forced circulation lead—bismuth cooled reactor. The results showed that a small 


amount of steam migrated into the core under the fracture condition of the lower head. 


Dinh’s study showed that when the steam generator is ruptured near the main pump, 
water vapor is bound to the core by lead—bismuth, and there is a large amount of 
water vapor in the active zoneľl. Such conditions cause changes in the reactivity of 
the fuel rods in the core, leading to the deterioration of heat transfer and endangering 
the safety of the reactor. Therefore, it is important to study the two-phase flow in fuel 
assemblies during SGTR accidents. Because it is difficult to use X-ray devices to 
study the internal flow field characteristics of fuel assemblies in lead-based reactors!?!], 
certain thermal—hydraulic parameters were obtained from numerical simulations. Gu 
et al.??] developed the multi-physics field coupling program MPC-LBE, which can 
perform thermal hydraulic analyses for LBE cooling-pool reactor types. Liu et al.!?! 
performed calculations for five turbulent Prandtl number models and accurately 
predicted the heat transfer characteristics of low-Prandtl number fluids in fuel 
assembly rod bundle subchannels. Yunxiang et al.!*4! used the k-æ SST model to study 
single-phase flow fields in lead-based reactors. Zhou et al.'*>! calculated the friction 
pressure drop inside a wire-wound rod bundle using eight different models, and the 
results showed that the friction coefficient was related to parameters such as the 
number of rod bundles and knot diameter ratio. Suzuki et al.'*°! studied the two-phase 
bubble flow in a lead—bismuth alloy. Wang et al.?7! developed a thermal hydraulic 
analysis program for lead—bismuth fast reactor ring fuel based on a closed parallel 
multichannel model and a quantum genetic algorithm. Lanting et al.'°*! numerically 
calculated the bubble flow in liquid heavy metals in a tube and analyzed the 
morphological changes in the bubbles as they rose in the tube. Jeltsov et al." studied 
the accumulation and flow of bubbles in a core during SGTR accidents, and the 
results indicated that small bubbles were more likely to remain in the core. 

In the development of lead-based fast reactors, SGTR accidents are important for 
analyzing severe accidents. The two-phase flow of the fuel assembly during an SGTR 
accident is directly related to the safety of the lead-based fast reactor design. Currently, 
there are few studies on two-phase flows in the fuel assemblies of lead-based reactors, 
and the experimental method is not ready for implementation. Therefore, the 


open-source OpenFOAM software was used to predict the two-phase flow of the fuel 
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assembly under SGTR accident conditions. This calculation can provide technical 
support for subsequent experimental methods and validate the design and 


development of lead-based fast reactors. 


2 Models and Methodologies 


2.1 Governing Equations of the LBEsteamEulerFoam Solver 

The Euler method assumes that each phase is on a mutually permeable 
continuum. Thus, each phase has independent physical properties such as the pressure, 
temperature, and velocity. It is necessary to analyze the interphase force to seal the 
governing equation when simulating a gas-liquid two-phase flow using the Euler 
method. In a liquid metal gas-liquid two-phase flow, the interphase mass transfer 
between the gas and liquid phases is generally ignored, and the governing equation is 


<(ap,)+¥-(apu)=0 at 


2 Capu) =-V-(a,puu)-a,Vp-V- |, (Vu +Vu ") a.p,g+F, (2) 
where ai, pi, and u; represent the macroscopic volume fraction, density, and velocity of 
phase i, respectively, and Fy represents the interphase force of the gas-liquid phase. Fy 
is primarily used to represent the transfer of momentum between phases, which is the 


result of the combined action of various forces”!: 


F, = F, +F, + Fy, + Fpp + Fyy (3) 


where Fp represents the drag force caused by the relative motions of the gas and 
liquid phases. FL represents the lift and lateral forces caused by the pressure 
difference perpendicular to the direction of bubble motion. Fw. stands for wall 
lubrication, which is generated by the slip velocity near the wall. Frp represents the 
turbulent diffusion force, which is the traction force caused by the liquid turbulence 
on the bubble. Fvm represents the virtual mass force, which is generated by changes in 
the relative acceleration. 
The drag force on a single bubble can be expressed as!*"!: 


u, u, 4; (4) 
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The shear-induced lift force of the bubble isP: 


F, =-C,p,a,u, x (V xu, ) (5) 


The mathematical expression of the wall lubrication force isP”: 


a 
Fy = Cy, A w 
d, (6) 
The turbulence dispersion force model is); 
Fp =- 4 ( _ -]va 
Pr, \ l-a a (7) 
The virtual mass force can be expressed as!*“); 
du, du 
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where Cp, p1, ur, and Aa represent the drag coefficient of the bubble, liquid density, 
relative velocity of the two phases, and projected area of the gas phase, respectively. 
Cı represents the lift coefficient. ag represents the volume fraction of the gas phase, 
and u represents the velocity of the liquid phase. Cwi is the wall lubrication force 
coefficient. Uw is the tangential component of the relative velocity at the wall. mw is 
the unit normal vector of the wall. Pra represents the turbulent Prandtl number with a 
discrete volume fraction, and yı represents the viscosity of the moving vortex in the 
liquid phase. Cvm is the virtual mass—force coefficient. 

The standard k-e model!*! is the most simple and economical method for 
modeling turbulence; it has good convergence and can accurately predict the flow in 
the pipeline, shear flow, etc. The turbulent kinetic energy k and turbulent kinetic 
energy dissipation rate ¢ are expressed as follows: 
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Then, Equations (9) and (10) are coupled by the following equations: 


2 
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(10) 


G, = wVU-| VU+(VU)’ | (12) 


where fm is the dynamic viscosity of a laminar fluid; yz; is the dynamic viscosity of a 
turbulent fluid; ox and o: represent the diffusion Prandtl numbers; the model constants 
Ca, Cer, and C2 are 0.09, 1.44, and 1.92, respectively. 

The Prandtl number of liquid lead-bismuth is nonlinear. The turbulence effect 
was considered in subsequent calculations; therefore, the Prandtl number of liquid 


lead—bismuth was modified? as: 


4.12, P, < 1000 (13) 
= OIP 
i= 0.01. __ 1000 < P < 6000 
[0.018P -7 -4| 
4.5, P, <1000 
A=45.4-9x107P,,1000< P, < 2000 
3.6, P, > 2000 


OpenFOAM provides users with many operators that can be used to solve partial 
differential equations. In this study, each physical quantity in the governing equation 
was first defined on the open-source platform, so that each quantity could be called 
for numerical calculation in the operation. Next, the governing equations of the 
physical quantities were combined using the OpenFOAM operators. OpenFOAM 
provides methods for explicit and implicit calculations that can be invoked in the 
conservation equations. The second-order discretization method was used in 
OpenFOAM to discretize the physical quantities and ensure calculation accuracy. The 
improved governing equation was encapsulated in OpenFOAM as a self-developed 
solver called LBEsteamEulerFoam. This solver can be used for liquid 
lead/lead—bismuth and water vapor two-phase flows. 

2.2 Meshing the CiADS Assembly 

In advanced reactors of conventional design with a lead/lead—bismuth alloy as 
the coolant, the fuel rods are usually arranged in a triangular manner to form 
hexagonal assemblies. The subchannels of the assemblies are divided into three types: 
internal, edge, and corner. In this study, the CiADS subcritical reactor fuel assembly 


parameters!” were used for the simulation. The simulation area was the active area of 
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the fuel assembly. The axial height of the studied section was 720 mm. There were 60 
fuel rods in the regular hexagonal fuel assembly, and one stainless-steel rod was 
located in the center. A stainless-steel rod used as a locking device was designed to 
prevent the assembly from floating in a liquid lead—bismuth alloy. The fuel rod 
diameter was 13.1 mm, and the pin pitch was 15.1 mm. 

Considering the accuracy of the two-phase flow calculations, a structured grid 
was used to discretize the model. The grid partitioning results are shown in Fig. la, 
and the boundary layer settings on the outer surfaces of the fuel rods are shown in Fig. 
1b. In OpenFOAM, the mesh was checked with the "checkMesh" command, and the 
results show that the mesh met the calculation requirements. The total number of cells 
was 2303680, and the maximum skewness was 0.544; the mesh orthogonal quality 


exceeded 0.7, and the Y-plus value was approximately 30. 


Fig. 1 Grid delineation and sensitivity analysis. a Grid division of the assembly; b Grid division 
of the boundary zones. 
2.3 Boundary Conditions Setting and Mesh-independent Analysis 
In numerical simulations, appropriate conditions must be set for each boundary 
according to the actual working conditions. The CiADS method [57-38] was used to set 


the boundary conditions listed in Table 1. A transient mode was adopted to simulate 
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the bubble behavior in the fuel assembly under different heating boundary conditions; 
the simulations were divided into three cases, as listed in Table 2. Currently, there are 
no reliable experimental data on the two-phase flow of lead-based reactor fuel 
assemblies. Therefore, the gas-liquid phase was assumed to be evenly mixed when the 
bubbles were transmitted through the main pump. The assembly inlet was set for 
uniform mixing of both the gas and liquid phases. The gas content is the fraction of 
the phase occupied by gas in the grid. The gas content rate was varied from 1% to 


30% to verify the impact of the bubble distribution, the velocity, and other operational 


parameters. 
Table 1 Boundary conditions for the CiADS fuel assembly 
Boundary Boundary conditions 
Average temperature of the inlet was set at 553.15K. 
Inlet 
The coolant velocity of the inlet was set at 0.36m/s. 
Outlet Pressure outlet with a gauge pressure of 0 Pab”?! 
60 fuel rods No-slip surface 
One stainless steel rod Adiabatic and no-slip surface 
Assembly box Adiabatic and no-slip surface 
Table 2 The setting of different cases 
Case number Heating boundary condition Coolant 
Case 1 Fixed heat flux 104483W/m? Liquid LBE 
Case 2 Fixed temperature 700K Liquid LBE 


Heat flux with cosine distribution in 
Case 3 Liquid LBE 
axial direction 24:37] 


To eliminate the influence of grids on the accuracy of the calculation results 
when solving partial differential equations, a grid-independence analysis should be 
carried out. Five groups of grids with different degrees of density were divided 
according to the fuel assembly model. The maximum temperature and velocity at the 
fuel assembly outlet were selected for the grid independence analysis. Fig. 2 shows 


the calculation results for the five groups of grids. The case with 2303680 grids 
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yielded similar results to that with 2786670 grids, while there was a large deviation 
for 1678372 grids. The calculated results for 5467372 and 12072264 grids are similar 
to those for 2303680 grids. However, the calculation times for 5467372 and 12072264 
grids were significantly longer. Considering the accuracy of the results and the 
calculation efficiency, the case with 2303680 grids was selected for subsequent 


calculations. 
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Fig. 2 Grid convergence study 

2.4 Coolant Properties 

A liquid lead—bismuth alloy was used as the coolant for the CiADS fuel assembly. 
Owing to the special thermal properties of Pb-Bi, the calculation process differed 
from that of conventional coolants. Therefore, the thermophysical parameters were 
experimentally measured and fitted to the correlations!*”! listed in Table 3. The steam 
properties obtained from the international standard IAPWS-IF97 steam table!*!! are 
defined in the OpenFOAM material database. 


Table 3 Properties of liquid LBE 


Properties Liquid LBE 
e[kg/m3] 11096 —1.3236T 
cpl J/(kg -K )] 159.0—0.0272T +7.12-10° (T)? 


u[Pa-s] 4.94 -10% -exp(754.1/T) 


27 


28 


A[W/(m-K )] 3.614+1.517-10°T -1.741-10° (TY 


2.5 Simulation Model Verification 

Owing to the opacity of liquid lead—bismuth and the difficulty in X-ray detection, 
experimental data on liquid lead—bismuth two-phase flows are scarce. To verify the 
accuracy of LBEsteamEulerFoam in the simulation of lead-based fast reactors, the 
single-phase flow of a liquid lead—bismuth cooled assembly containing spacer wires 
was simulated. The LBEsteamEulerFoam solver was verified by setting the inlet gas 
content to 1x10 without closing the inter-phase force factor. This order of magnitude 
is attributed to the oxygen control level of the LBE loop. In 2016, KIT conducted a 
19-rod bundle wire-wrapped positioning fuel assembly heat transfer experiment using 
LBE as a coolant. The experimental setup was installed in the vertical section of the 
THEADES loop. The specific boundary conditions used in the experiment were as 
follows: inlet temperature Tin=473 K; inlet mass flow rate M=19.18 kg/s; and total 
heating power Q=197.04 kW. The heating power was evenly distributed over the 
heating section of each rod. According to the MYRRHA reactor design, the 
experimental data from Pacio!*?! and the calculated results of the experimentally 
checked SACOS-PB subchannel program!!! were compared under the same 
conditions. The SACOS-PB subchannel program is a necessary analysis tool for 
reactor thermal-hydraulic design and safety analysis. To enhance the calculation 
efficiency, the inlet and outlet boundaries were set as periodic flow boundaries 
according to the periodicity of the spacer wires. Pacio selected three measurement 
positions with different axial heights in the heating section. Because the gas-liquid 
two-phase flow was more affected by the change in the axial height, the experimental 
data at an axial height of 820 mm were selected for comparison. The measurement 
points were arranged in five subchannels: 3, 14, 15, 29, and 39. The calculated results 
for the center positions of these five subchannels were extracted for comparison with 
the experimental results. A comparison of the coolant temperatures calculated for 
different subchannels at the axial height of the active zone of 820 mm is shown in Fig. 
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Fig. 3 Coolant temperature distribution at an axial height of 820 mm in the active zone. 

As shown in Fig. 3, the coolant temperatures in different subchannels on the 
same cross section are different. The coolant temperature in the inner channels near 
the central rod bundle was higher. The inner channel coolant temperature is higher 
than the edge and corner channel coolant temperatures. The coolant temperatures of 
the SACOS-PB and LBEsteamEulerFoam solvers in different subchannels in the same 
section were in good agreement under the same calculation conditions, and the 
maximum relative error was less than 1.13%. Compared to the experimental data, the 
maximum relative error of the LBEsteamEulerFoam solver was less than 4.17%. 
Therefore, the simulation results of the LBEsteamEulerFoam solver are accurate, and 
the calculation accuracy is high for 10° gas content. Hence, this solver can be used 


for the subsequent simulation analysis of two-phase flows. 


3 Results and discussion 


3.1 Bubble Gathering 

Fig. 4 shows the bubble distribution scenarios in Case 1 when the water vapor 
content at the fuel assembly inlet was 10%. Owing to the transient nature of the 
multiphase flow calculation, three typical bubble distribution scenarios were selected 
for analysis in the middle and outlet sections. Data from the three time points were 
extracted. Three sets of data from the same moment were normalized. The error bars 


represent fluctuations in the calculated data. As shown in Fig. 4, the bubbles 
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accumulated at the edge and corner channels of the fuel assembly. The fraction of 
bubbles in the edge and corner channels increased with time. The distribution results 
provide a reference for the aggregation analysis of liquid heavy metals in fuel 
assemblies. There is a large temperature distortion in the bubble-gathering area; 
therefore, it is necessary to further improve the high-temperature resistance of the fuel 
assembly box to cope with a SGTR accident. The degree of bubble aggregation in the 
edge channel of the outlet exceeded that in the middle section, which indicates that 
the bubbles aggregated when rising in the assembly. 

Because the corner channel consists of two walls at 120°, the flow in this region 
is influenced by both the walls of the fuel assembly cartridge and fuel rods. The 
heating area in this region was smaller than those in the other channels, complicating 
the flow characteristics in this region and resulting in uneven coolant velocity 
distributions in the same cross section. The coolant velocity was reported to be lower 
in the corner channels. As the coolant flowed through the fuel assembly, the bubbles 
were dragged by the coolant and gradually deflected from the high-velocity region to 
the low-velocity region. Therefore, bubbles accumulate at the edge and corner 


channels of the fuel assembly. 


Fig. 4 Bubble distribution in Case 1 with 10% inlet steam content. 


3.2 Outlet Section Velocity 

Fig. 5 shows the maximum velocity of the two-phase fluid in the fuel assembly 
outlet section. Based on the three different cases, 30 sets of data with gas contents 
ranging from 1% to 30% were calculated. Because of the low density of water vapor 
and its expansion effect, the velocity of the two-phase flow is closely related to the 
imported water vapor content. The figure shows that the maximum velocity of the 
liquid metal at the fuel assembly outlet increased with the inlet water vapor content. 
When the inlet water vapor content was less than 15%, the fluid velocity fluctuation at 
the assembly outlet was small, and the growth rate was slow. When the inlet water 
vapor content exceeded 15%, the two-phase flow velocity fluctuated significantly, and 
the fluctuation range increased gradually. This phenomenon indicates that the 
presence and migration of the gas phase disturbed the coolant flow in the fuel 
assembly. A comparison with Fig. 5 shows that the heating condition with a fixed 
temperature causes the component outlet to obtain a larger two-phase flow velocity 
under a high gas content. When the axial power distribution was set to a cosine shape, 
larger velocity fluctuations were observed at the outlet. 

Owing to the expansion of steam bubbles, the velocity of the two-phase flow is 
related to the inlet steam content. The gas phase is less intensive, and large bubbles 
form along the axial direction, occupying a larger cross-sectional area. Because 
lead—bismuth and steam bubbles may have no phase conversion, the cross-sectional 
area of the liquid phase may decrease. Owing to the continuity principle, the 
liquid-phase velocity may increase. The velocities of these two phases increase by 
different magnitudes. The Euler method employs two sets of NS equations and can 
solve for both phases separately; thus, it is more suitable for simulations of 
high-gas-content two-phase flows. The lift force between the two phases may increase 
the velocity of the gas phase. However, the growth of the gas-phase velocity is 


smaller than that of the liquid phase because of the larger resistance. 


20 


21 


—#— Fixed heat flux value 
—e— Fixed temperature 
—A— Cosine distribution 


Maximum Coolant Velocity/(m/s) 


0 5 10 15 20 25 30 


Steam volume fraction (%) 
Fig. 5 Maximum coolant velocity at the assembly outlet section. 
3.3 Temperature Distributions at the Middle and Outlet Sections 

During an SGTR accident, the coolant temperature in the fuel assembly is 
affected by the water vapor content. The design parameters of the fuel assembly, such 
as the material selection or geometrical design, may significantly depend on the 
anticipated operating temperature. Fig. 6a shows the maximum coolant temperature 
change in the axial intermediate section of the fuel assembly when the inlet water 
vapor content changes from 1% to 30%. Because the temperature change of the 
coolant is not obvious for fixed-temperature heating, heating conditions with a fixed 
heat flux and cosine distribution were selected for the analysis. When the inlet water 
vapor content increased from 1% to 15%, the maximum coolant temperature in the 
middle section increased slowly, and the water vapor content did not dominate the 
two-phase flow. When the imported water vapor content exceeded 15%, the coolant 
temperature increased significantly and obvious fluctuations occurred, which gives 
insight into the degree of heat resistance fatigue of the fuel rods. 

Fig. 6b shows the relationship between the inlet vapor content and temperature at 
the outlet section with fixed and cosine-shaped heat fluxes. When the inlet water 
vapor content was less than 15%, the coolant temperature in the assembly remained 
within the designed safety range of the fuel assembly. The coolant temperature 


fluctuation was small, indicating that the fuel assembly could tolerate small amounts 


of water vapor under normal operating conditions. When the vapor content at the inlet 
exceeded 15%, the maximum temperature of the coolant at the outlet section 
increased noticeably and the fluctuation became more significant. Thus, the fuel 
assembly can tolerate greater thermal stress and fatigue. 

Fig. 6c illustrates the temperature distribution of fuel rods 2 and 40 along the 
axial direction under heating with a cosine-distributed heat flow density. The 
difference between the single-phase flow temperatures of fuel rods 2 and 40 was 
minimal because there was no cross mixing of the spacer wires. In the two-phase flow 
of the fuel assembly, the axial temperatures of both fuel rods increased with the axial 
height. However, the maximum temperature increase on the cladding surface of fuel 
rod 2 exceeded that of fuel rod 40. The largest temperature growth rate of the cladding 
surface was observed at Z=0.36 m, where the line power density reached its 
maximum. The temperature at the outlet was close to 850 K for fuel rod 2, which was 
considerably higher than the maximum temperature obtained by the CiADS 
calculations®7!, This calculation was instantaneous, and the location of the grid with 
the maximum gas content changed as the calculation time increased. The contour 
results demonstrate that the bubbles gathered in the edge and corner channels. The 
hottest locations of the assembly outlet were typically distributed in the inner 
channels. 

Owing to the relatively low specific heat capacity of steam, the cooling capacity 
of the cladding in the gas phase is much lower than that of liquid lead—bismuth, which 
may cause unexpected temperature increases in the fuel assembly. Simultaneously, 
owing to the narrow subchannel area within the fuel assembly, large bubbles may 
block the subchannel near the assembly outlet, causing heat transfer deterioration. 
When the heat flow density was distributed along the axial cosine, the maximum 
linear power rating was observed at Z=0.36 m. Because bubbles may accumulate at 
this location, the heat transfer coefficient between the cladding surface and coolant 
may deteriorate significantly. Because of the relatively slim cladding wall, the heat 


capacity was limited, which may have caused partial overheating of the cladding. 
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Fig. 6 Maximum temperatures. a Maximum LBE temperature at the assembly middle section; 


b maximum LBE temperature at the assembly outlet section; ¢ Comparison of the maximum 


cladding temperature distribution along the axial direction of the fuel rods. 
3.4 Maximum Steam Volumetric Fraction 

Fig. 7a—c shows the relationship between the maximum vapor content in the 
middle and outlet sections of the fuel assembly and the vapor content in the inlet 
section. The maximum gas holdup of the grid at the outlet section of the fuel assembly 
exceeded that in the middle section. According to the contours in Section 3.1, the 
degree of bubble accumulation in the edge channels at the fuel assembly outlet 
exceeded that at the edge channels in the middle section. These results indicate that 
the bubbles accumulated during rising in the fuel assembly subchannels. When the 
gas content of the inlet was less than 2%, the maximum gas contents of the grid at the 
outlet and middle sections of the fuel assembly were slightly different. This indicates 
that the bubble accumulation phenomenon was not obvious at this time, and no large 
bubbles were generated. When the vapor content at the inlet increased from 5% to 
15%, the grid gas holdup at the outlet section increased steeply, indicating that the 
maximum vapor content of the grid increased more rapidly during this process. 
Moreover, the gas-holdup growth slope of the grid in the outlet section was higher 
than that of the grid in the middle section, indicating that large bubbles gradually 
accumulated in the upper part of the fuel assembly. When the inlet gas content 
exceeded 15%, the growth rate of the maximum gas content in the grid was modest, 
and the maximum gas content in the grid gradually reached its maximum value. When 
the gas content at the inlet exceeded 25%, there was little difference between the 
maximum gas content at the outlet and middle sections of the fuel assembly. This 
indicates that the bubble coalescence process in the upper part of the assembly 
became saturated at a higher gas content. The percentage of gas content in the grid 
increased owing to the presence of large air bubbles, and the maximum gas content of 


the grid increased with increasing axial height. 
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Fig. 7 Maximum steam content of the grid. a Case 1; b Case 2; ¢ Case 3. 
3.5 Comparison with the VOF model 

At present, there are few studies on two-phase flows in the fuel assemblies of 
lead-based fast reactors, and the experimental method is not yet ready for 
implementation. Multiphase flow simulations are divided into an interface class 
model (VOF method) and a high-phase-fraction multiphase flow model. The VOF 
model is a surface-tracing method that can be used to determine the interface between 
a variety of mutually incompatible fluids. Sheng et al.'*! used the VOF method to 
simulate the two-phase flow of an ADS fuel assembly with different inlet water vapor 
contents, heating boundary conditions, and coolants. The results show that water 
vapor accumulated heat in the fuel assembly, and the temperature fluctuated 
significantly. Because of the lack of experimental data, the Euler calculation results 
were compared with the VOF calculation results. Fig. 8a shows the accumulation of 
water vapor bubbles at the edge and corner channels calculated using the Euler and 
VOF methods. The gas bubbles in the corner channel of the fuel assembly 
accumulated into clumps. The degree of bubble aggregation in the corner channel 
exceeded that in the edge channel. According to the contour comparison, larger 
bubbles were formed in the corner channel using the VOF method, and the degree of 
bubble aggregation in the side channel exceeded that calculated using the Euler 
method. The degree of bubble aggregation significantly influenced the coolant speed 
and temperature. 

Fig. 8b shows the maximum temperature of the fuel assembly outlet calculated 
using the Euler and VOF methods. The maximum temperature of the component 
outlet increased with the gradual increase in the inlet water vapor content. When the 
imported vapor content exceeded 15%, the coolant temperature increased significantly. 
The results of the two solvers agree well with the upward trend, and the maximum 
relative error of the temperature was less than 5.83%. The simulation time of the 
LBEsteamFoam solver 10 s after the accident was 42.3 h, whereas that of the 
LBEsteamEulerFoam solver was only 1.64 h with the same hardware resources. 


Compared with the VOF method, the Euler method has a higher calculation efficiency 
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when dealing with two-phase flows in LFR fuel assemblies. There is a significant 
difference between the two methods in terms of the number of meshes required. The 
VOF method is a direct simulation of a multiphase flow, similar to the DNS. Because 
the VOF method must capture the boundary between the two phases, a finer mesh 
may be required. The Euler method, however, need not finely grasp the boundary 
between the two phases and therefore requires a smaller mesh. Thus, a higher 
calculation efficiency can be expected using the Euler method. 

The VOF method can capture the boundary between two phases and provide a 
good description of bubble accumulation in the corner and edge channels of the fuel 
assembly. Each phase in the Euler method has its own physical properties, which can 
better describe parameters such as the bubble diameter, interphase force model, and 
multi-scale distribution. The Euler method has a unique advantage for 
high-phase-fraction multiphase flow calculations. Owing to the deviation in the 
calculation results of the Euler and VOF methods, experiments are required for 
further verification. The fuel assembly test equipment requires assembly boxes 
because bubbles accumulate at the edge and corner channels. Simultaneously, 
measurement points were arranged in the corner and edge regions near the assembly 


box to monitor the water vapor phase in the two-phase flow. 
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Fig. 8 Comparison with the VOF method. a Steam bubbles at the periphery regions; b 


Comparison of the maximum coolant temperatures at the outlet. 


4 Conclusions 


Improving the safety of lead-based fast reactors is beneficial when studying 
pipeline rupture accidents in steam generators. In this study, the open-source 
computational fluid mechanics software OpenFOAM was used to calculate the 
two-phase flow of liquid heavy metal and water vapor in the fuel assembly. By 
changing heating boundary conditions, the bubble aggregation, velocity, temperature 
and maximum gas content of the cross-sectional grid in the fuel assembly under 
different working conditions were analyzed, and the conclusions were obtained as 
follows: 

1) At the fuel assembly outlet, bubbles are more likely to accumulate in the edge 
and corner channels near the assembly box, and the degree of accumulation in the 
corner channels is greater than that in the edge channels. As the inlet water vapor 
content increases, the two-phase flow velocity at the fuel assembly outlet may 
increase. 

2) As the inlet vapor content increases, the temperature inside the fuel assembly 
will gradually increase. Due to the weak heat-conducting capability of water vapor, 


the maximum coolant temperature that can be reached in the fuel assembly might 


increase. When the inlet water vapor content is less than 15%, the maximum 
temperature at the fuel assembly outlet is still within the design range. When the inlet 
water vapor content is more than 15%, the maximum temperature at the fuel assembly 
outlet rises sharply, which will easily reach the working limit of fuel rods and thus 
affect the safety behavior of fuel assemblies. At the same time, when the inlet vapor 
content is set higher than 15%, the temperature fluctuation may become larger, which 
will aggravate the thermal fatigue stress of the fuel rod and the assembly box. 

3) With the increase of inlet steam content, the maximum gas content of the grid 
at the fuel assembly outlet also increases gradually. When the inlet water vapor 
content is less than 5%, the gas phase in the fuel assembly is not dominant. When the 
inlet water vapor content is in the range of 5% to 20%, the upper part of the fuel 
assembly will gradually accumulate to form large bubbles. A large number of bubbles 
coalesce and break, so the maximum gas content of the grid at the fuel assembly 
section increases significantly. When the gas content of the inlet is greater than 20%, 
the maximum gas content of the grid no longer shows an obvious upward trend. 

4) Compared with the VOF method, the Euler method has higher calculation 
efficiency when dealing with the two-phase flow in LFR fuel assemblies. In addition, 
each phase in Euler method has its own physical properties, which can be considered 
to provide better descriptions of parameters such as bubble diameter, interphase force 
model, and multi-scale distribution. However, the Euler method may cause the 
underestimation of the void fraction and lead to the underestimation of the positive 
reactivity feedback of the lead-based reactors, so it is still necessary to cross check 


with two-phase flow experimental data in the future. 
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